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ABSTRACT 

We use numerical simulations of turbulent cluster-forming regions to study the nature of dense 
filamentary structures in star formation. Using four hydrodynamic and magnetohydro dynamic sim¬ 
ulations chosen to match observations, we identify filaments in the resulting column density maps 
and analyze their properties. We calculate the radial column density profiles of the filaments every 
0.05 Myr and fit the profiles with the modified isothermal and pressure confined isothermal cylinder 
models, finding reasonable fits for either model. The filaments formed in the simulations have similar 
radial column density profiles to those observed. Magnetic fields provide additional pressure support 
to the filaments, making ‘puffier’ filaments less prone to fragmentation than in the pure hydrodynamic 
case, which continue to condense at a slower rate. In the higher density simulations, the filaments 
grow faster through the increased importance of gravity. Not all of the filaments identified in the 
simulations will evolve to form stars: some expand and disperse. Given these different filament evo¬ 
lutionary paths, the trends in bulk filament width as a function of time, magnetic field strength, or 
density, are weak, and all cases are reasonably consistent with the finding of a constant filament width 
in different star-forming regions. In the simulations, the mean FWHM lies between 0.06 and 0.26 pc 
for all times and initial conditions, with most lying between 0.1 to 0.15 pc; the range in FWHMs 
are, however, larger than seen in typical Herschel analyses. Finally, the filaments display a wealth of 
substructure similar to the recent discovery of filament bundles in Taurus. 

Subject headings: 


1. INTRODUCTION 

Filaments appear to be an important ingredient in the 
formation of stars. While filaments have been known to 
be associated with star forming regions for decades (e.g., 
Schneider & Elmegreen 1979; Bally et al. 1987), observa¬ 
tions from the Herschel Space Telescope^ particularly the 
Gould Belt (Andre et al. 2010) and HOBYS (Motte et al. 
2010) Legacy Surveys have underlined the prevalence of 
filamentary structures within star forming regions. With 
HerscheFs unprecedented ability to sensitively map large 
areas of the sky, several common properties of filaments 
have now been identified. First, filaments appear to not 
be well represented by the Ostriker (1964) equilibrium 
model of an isothermal cylinder; the column density pro¬ 
file is shallower (e.g., Arzoumanian et al. 2011). This may 
indicate that magnetic fields (e.g., Fiege & Pudritz 2000) 
contribute to supporting the filament from collapse, al¬ 
though Smith et al. (2014) demonstrate that filaments 
formed in purely turbulent environments also have a sim¬ 
ilarly shallow slope. Rotation may also lead to a shal¬ 
lower slope (Recchi et al. 2014). Second, the mass per 
unit length of filaments appears to correlate with star- 
formation activity: filaments with mass per unit length 
less than the value needed for collapse of an isothermal 
cylinder (Ostriker 1964) tend to be associated with re- 
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gions which are forming few if any stars, while filaments 
with supercritical mass per unit length values tend to be 
associated with active star forming regions (e.g., Arzou¬ 
manian et al. 2011; Hennemann et al. 2012). What is still 
unclear, however, is what forces dominate the formation 
and evolution of the filaments, and how the filaments 
contribute to star formation. For example, are the fil¬ 
aments formed primarily through turbulent shocks, or 
under the influence of magnetic fields or gravity? Does 
turbulence control the ability of filaments to fragment 
into star-forming cores? What forces set the observed 
(column) density profiles? And do filaments primarily 
provide a denser collection of gas to promote local star 
formation (e.g. Hacar & Tafalla 2011), or do they play 
a significant role in providing a conduit of mass for the 
formation of larger stellar clusters, which appear to form 
preferentially at the intersection of several filaments (see 
e.g., Myers 2009, 2011; Schneider et al. 2012; Hennemann 
et al. 2012; Kirk et al. 2013)? 

In this paper, we investigate the first of these issues, 
namely the formation and evolution of filaments, through 
the analysis of our numerical simulations. We compare 
the column density properties of filaments formed within 
four different simulations: higher and lower density, and 
with and without magnetic fields. These analyses provide 
a complementary look at simulations to those recently 
published in Smith et al. (2014), where the influence 
of different types of turbulence on filament properties 
was examined, but the effect of the inclusion of magnetic 
fields or differing initial mean densities was not. 

We find that while the largest-scale structures in the 
gas are set by turbulent motions, and appear similar in 
all four simulations, magnetic fields and gravity do in¬ 
fluence the properties of individual filaments. In par- 
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ticular, magnetic fields cushion the initial turbulent gas 
compressions, leading to filaments which are initially less 
condensed, and subsequently evolve more slowly (due to 
the weaker gravitational pull) than the corresponding hy¬ 
drodynamic case. We note that the simulations we an¬ 
alyze were only able to be run for a few tenths of the 
global free-fall time, limiting our sensitivity to later-time 
evolutionary trends. The simulated filaments have prop¬ 
erties which are consistent with those measured in real 
filaments characterized by Herschel^ suggesting that the 
general insights gained with these simulations are appli¬ 
cable to real molecular clouds. Finally, turbulence and 
magnetic fields, and not just the thermal properties of 
molecular gas, appears to set the critical conditions for 
gravitational instability leading to star formation. 

In what follows, we first discuss our numerical meth¬ 
ods and simulations (Section 2), discuss the basic fila¬ 
ment properties resulting from the simulations (Section 
3), compare various models of filament structure (Sec¬ 
tion 4), and examine the effects of spatial resolution in 
characterizing filaments (Section 5). We discuss our re¬ 
sults and their implications, as well as the limitations of 
our present analysis in Section 6. 

2. NUMERICAL METHODS 
2.1. Simulation Setup 

We used the flash hydrodynamics code (Fryxell et al. 
2000) version 2.5 to perform numerical simulations of 
molecular clumps, i.e., parsec-scale condensations of gas 
capable of forming a cluster of stars, flash solves 
the fluid-dynamical equations on an adaptive Eulerian 
grid, making use of the paramesh library (Olson et al. 
1999; MacNeice et al. 2000). It includes self-gravity, La- 
grangian sink particles to represent gravitationally col¬ 
lapsing cores and (proto)stars (Banerjee et al. 2009; Fed- 
errath et al. 2010), and gas cooling by dust and by molec¬ 
ular lines (Banerjee et al. 2006). Stellar properties are 
self-consistently evolved via a one-zone model (Offner 
et al. 2009; Klassen et al. 2012). 

We initialize our simulation volume with a turbulent 
velocity field. The turbulence is a mixture of compres¬ 
sive and solenoidal turbulence (Federrath et al. 2008; 
Girichidis et al. 2011) with a Burgers spectrum, i.e. 
Ek oc k~‘^ as in Girichidis et al. (2011), and largest modes 
having a size scale roughly equal to the side length of the 
simulation box. See also Larson (1981); Boldyrev (2002); 
Heyer & Brunt (2004). The turbulent velocity field has 
a root-mean-square Mach number of 6. 

We perform a grid of simulations in a cube-shaped vol¬ 
ume containing either approximately 500 Mq or 2000 Mq 
of molecular gas with a power-law density profile scaling 
as p{r) = pcr~^^^. The choice of density profile is moti¬ 
vated by observations of dense gas associated with high- 
mass star formation (Pirogov 2009); Kauffmann et al. 
(2010) similarly analyze a suite of dust emission and ex¬ 
tinction maps of molecular clouds within the solar neigh¬ 
bourhood, and find that those which are not forming 
high-mass stars obey p{r) oc The simulation vol¬ 

ume has a side length of 2 pc, and the molecular gas is 
at an initial temperature of 10 K. 

These initial conditions were chosen to be represen¬ 
tative of nearby molecular clumps, with a focus on 
NGG1333, a cluster-forming region within the Perseus 


molecular cloud, located roughly 250 pc away, and cur¬ 
rently forming a young cluster of low- and intermediate- 
mass stars (Walawender et al. 2008). Using a large-scale 
column density map derived from 2MASS-based extinc¬ 
tion, Kirk et al. (2006) estimate that NGG1333 contains 
^1000 Mq within a radius of ^1 pc; the simulations con¬ 
tain 500 and 2000 Mq within a 2 pc cube, thus bracketing 
NGG1333’s mean density. The free-fall time for these 
simulations is ^ 1 and 0.5 Myr respectively. A Mach 
number of 6 is consistent with the typical ^^GO velocity 
dispersion measured across NGG1333 reported in Kirk 
et al. (2010), and we also note is also consistent with 
the standard linewidth-size relationship Larson (1981). 
Molecular clumps tend to have temperatures of 10-20 K 
(Bergin & Tafalla 2007), and pointed observations to¬ 
ward dense cores in Perseus (Rosolowsky et al. 2008) 
have a mean temperature of 11 K, although those found 
in NGG1333 and other clustered environments tend to 
have slightly higher values (Schnee et al. 2009; Foster 
et al. 2009). Similarly, the dust temperature is estimated 
to be slightly elevated in areas near luminous young pro¬ 
tostars (Hatchell et al. 2013). None of these heating ef¬ 
fects, however, would have been present prior to the onset 
of star formation in the region, suggesting that an initial 
temperature of 10 K is reasonable. 

We used the same initial turbulent velocity field for 
each simulation, but compared magnetohydrodynamic 
runs with pure hydro simulations where the magnetic 
field strength was set to zero. When including mag¬ 
netic effects, we initialize a magnetic field parallel to 
the z-axis with uniform field strength. We select a field 
strength for our MHD simulation so our mass-to-flux ra¬ 
tio is A ^ 1 — 2; this is slightly stronger than the typical 
range estimated by Grutcher et al. (2010) of 2 — 3. The 
mass-to-flux ratio is given by 


Mtot VG 

7ri?2 < B > 0.13 


( 1 ) 


where Mtot is the total cloud mass, R the cloud radius, 
and < B > the initial mean magnetic field strength. 
The factor of 0.13 is required to normalize the flux ratio 
relative to the critical value where the magnetic field just 
prevents gravitational collapse (Mouschovias & Spitzer 
1976; Seifried et al. 2011). High-mass star forming cores 
typically have values A < 5 (Falgarone et al. 2008; Girart 
et al. 2009; Beuther et al. 2010). 

Table 1 lists the parameters for the grid of simulations 
run. We note that while stars (sink cells) do form in 
all of our simulations, as we would expect in reality, the 
resolution (50 AU) is insufficient to correctly predict the 
masses of the stars that form; tests we ran with an in¬ 
creased resolution led to a larger number of lower mass 
stars. This is not a problem for our analysis, as the reso¬ 
lution is more than sufficient to characterize the structure 
of the filamentary gas at observable scales. Furthermore, 
the simulations are stopped at an early enough time that 
stellar feedback would not have had time to influence the 
evolution of the gas. 


2.2. Filament Identifieation 

The initial turbulent velocity field quickly results in a 
highly filamentary structure, as illustrated in Figure 1. 
We run each of our simulations until the filamentary 
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TABLE 1 

Simulation parameters 


Physical simulation parameters 

Parameter 



500HYD 

500MHD 

2000HYD 

2000MHD 

cloud radius 

[pc] 

Rq 

0.99978 

0.99978 

0.99978 

0.99978 

total cloud mass 

[M©] 

Aftot 

502.603 

502.603 

2152.11 

2152.11 

mean mass density 

[g/cm3] 

ip) 

4.256 X 10-21 4.256 x lO-^i 

1.822 X 10-20 

1.822 X 10-20 

mean number density 

[cm-3] 

{n) 

1188.98 

1188.98 

5091.14 

5091.14 

mean molecular weight 


P 

2.14 

2.14 

2.14 

2.14 

temperature 

[K] 

T 

10 

10 

10 

10 

sound speed 

[km/s] 

Cs 

0.196 

0.196 

0.196 

0.196 

rms Mach number 


M 

6.01 

6.01 

6.01 

6.01 

rms turbulent Alfvenic Mach Number 


Ma 

2.1 

2.1 

2.2 

2.2 

mean freefall time 

[Myr] 


0.74 

0.74 

0.370 

0.370 

sound crossing time 

[Myr] 

tsc 

9.96 

9.96 

9.96 

9.96 

turbulent crossing time 

[Myr] 

he 

1.66 

1.66 

1.66 

1.66 

Jeans length 

[pc] 

Aj 

0.413 

0.413 

0.199 

0.199 

Jeans volume 

[pc^] 

Vj 

0.294 

0.294 

0.033 

0.033 

Jeans mass 

[Mq] 

Mj 

4.42 

4.42 

2.13 

2.13 

magnetic field 

[mG] 

B 

0 

56.7 

0 

120.5 

mass-to-fiux ratio 


A 

oo 

1.17144 

oo 

2.35979 

rigid rotation angular frequency 

[rad/s] 

f^rot 

1.114e-14 

1.114e-14 

1.114e-14 

1.114e-14 

rotational energy fraction 


drot 

1.8 % 

1.8 % 

0.4% 

0.4% 

Numerical simulation parameters 

simulation box size 

[pc] 

hbox 

1.99956 

1.99956 

1.99956 

1.99956 

simulation box volume 

[pc"] 

Ebox 

7.99471 

7.99471 

7.99471 

7.99471 

smallest cell size 

[AU] 

Ax 

50.3465 

50.3465 

50.3465 

50.3465 

Simulation outcomes 

final simulation time 

[kyr] 

^final 

179.3 

232.2 

42.8 

49.6 

number of sink particles formed 


^sinks 

16 

6 

45 

3 

max sink mass 

[M©] 


2.01528 

9.53198 

19.4442 

31.2937 

min sink mass 

[M©] 


0.0264274 

0.164572 

0.00759937 

8.42947 

mean sink mass 

[M©] 


0.696485 

2.84511 

0.680372 

19.8616 

median sink mass 

[M©] 


0.525414 

1.56426 

0.0548092 

19.8616 


structure is well-developed; the simulation is stopped at 
0.2 to 0.3 free-fall times for the 500 Mq simulations (for 
the MHD and HD simulations respectively), and 0.13 
free-fall times for the 2000 Mq simulations. As flash is 
an adaptive mesh refinement (AMR) code, we first take 
the output files and map them to a uniform grid, down- 
sampling somewhat to allow the entire grid to fit into 
memory. Even with the downsampling, our resolution is 
^0.002 pc, much better than achievable with Herschel for 
nearby star-forming regions. We then project the den¬ 
sity along each of the coordinate axes to create column 
density maps. 

Figure 1 shows the column density in the X projec¬ 
tion for both the 500 Mq and 2000 Mq simulations at 
all time steps analyzed. Note that the MHD simulation 
was run for 0.15 Myr, while the HD simulation was run 
for 0.2 Myr for the 500 Mq simulations, giving one ad¬ 
ditional time step for our HD analyses. In this figure, all 
the panels have the same dynamic range shown for the 
greyscale column density, highlighting that material ac¬ 
cumulates into filamentary structures quite quickly (top 
and middle panel from left to right), and that having an 
initially higher density more rapidly leads to dense fila¬ 
mentary structures due to the increased importance of 
gravity (bottom row, left and middle panels). Finally, 
the presence of a magnetic field acts to slow the accu¬ 


mulation of material into dense filaments, as can be seen 
comparing the top and middle row panels, or the bottom 
row left and middle panels. We will return to this point 
in more detail in Section 3 and beyond. 

To extract the filamentary structure evident in Fig¬ 
ure 1 , we use the DisPerSE filament-finding algorithm^ 
described in Sousbie (2011) and Sousbie et al. (2011). 
The DisPerSE algorithm identifies persistent topolog¬ 
ical structures such as peaks, voids, and filaments, and 
is effective even if the image is noisy. It has been ex¬ 
tensively used on Herschel observations for filamentary 
structure identification, e.g. Arzoumanian et al. (2011); 
Schneider et al. (2012); Peretto et al. (2012); Palmeirim 
et al. (2013). In DisPerSE, there are several user- 
defined parameters to control the resulting filamentary 
network: persistence and robustness thresholds, smooth¬ 
ing, and a maximum angle. The two thresholds can be 
thought of as very roughly corresponding to criteria for a 
minimum absolute brightness (persistence threshold) and 
a minimum relative brightness compared to neighbour¬ 
ing features (robustness threshold). Smoothing removes 
small-scale ‘wiggles’ from the initial filament spine, while 
the angle is used to specify the minimum angular rota¬ 
tion between two initial filament spine segments that can 

® http://www2.iap.fr/users/sousbie/ 
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be joined together and still be classified as the same fil¬ 
ament. Filament spine segments which meet at a right 
angle, for example, are likely not part of the same fila¬ 
ment. 

We identify filaments using a persistence threshold of 
0.025 g cm“^ and a robustness threshold of 0.05 g cm“^ 
in the 500 Mq simulation (or 7 and 14 xlO^^cm”^) 
and thresholds of 0.1 g cm“^ and 0.2 g cm“^ (or 2.8 
and 5.6x10^^ cm“^) respectively for the 2000 Mq sim¬ 
ulation, smoothing the resulting filaments 1000 times, 
and allowing the initially identified filament segments to 
be connected for angles of less than 60 degrees (rela¬ 
tive to a straight line). These parameters were chosen 
after testing a range of values to determine which val¬ 
ues produced a filamentary structure that best matched 
visually-apparent structures. All of these thresholds for 
DisPerSE are above the standard ‘threshold for star 
formation’ found in nearby molecular clouds of around 
^ 5 —7x 10^^ cm“^ (e.g., Johnstone et al. 2004; Konyves 
et al. 2013). Unlike the Herschel analyses, we applied 
DisPerSE directly on the column density map. Since 
our column density maps include only the gas from the 
simulated star-forming clump, with no potential con¬ 
tribution from other dense structures within the larger 
cloud, we have less need than with Herschel data to apply 
filament-enhancing algorithms. Finally, we excluded sev¬ 
eral very short filaments that DisPerSE initially identi¬ 
fied - in order to accurately determine the filament profile 
(below), we set a minimum length of 0.1 pc. 

Figure 2 shows the network of filaments identified in 
the X projection of the 500 Mq and 2000 Mq HD simu¬ 
lations overlaid on their column density maps. 

One of the goals of our analysis is to track the time 
evolution of and the effect of magnetic fields on individ¬ 
ual filaments. In order to do so, DisPerSE was not 
used to identify a different network of filaments at every 
time step and magnetic field value, as this could poten¬ 
tially lead to different filaments being identified at dif¬ 
ferent snapshots. Instead, for each of the three projec¬ 
tions, we started with the network of filaments identified 
with DisPerSE at 0.15 Myr in the HD simulation, and 
then searched for the corresponding structures at dif¬ 
ferent times and with magnetic fields present. For the 
2000 Mq simulations, we instead started with the sin¬ 
gle 0.05 Myr time step. We started with an automated 
procedure to identify equivalent filaments at other time 
steps and / or with magnetic fields, by searching for local 
column density maxima near the reference set of filament 
spines. After this step, all filament spines were verified 
and adjusted as necessary by hand, using a combina¬ 
tion of visual inspection of the current column density 
snapshot and a movie of the time evolution of the col¬ 
umn density map for the 500 Mq simulations. The sim¬ 
ulations, particularly without the moderating presence 
of magnetic fields, form significant substructure on all 
scales, making it difficult to impossible for an automated 
procedure to correctly ‘follow’ the filaments in time and 
across initial conditions. 

There are several cases where a filament could not be 
fully traced to earlier times or in the corresponding sim¬ 
ulation with magnetic fields. Some, but not all, of these 
cases appear to be attributable to structures which are 
only apparent as filaments in 2D due to a coincidence of 


independent 3D structures; at other time steps, the real 
3D structures have moved by different amounts and no 
longer appear connected. We include these structures in 
our analysis where they do appear as a single filamen¬ 
tary structure, as any real observation which only has 
column density information is fallible to the same line of 
sight coincidence confusion. We will address the full 3D 
nature of filaments in these simulations in an upcoming 
paper. 

A comparison of Figure I and Figure 2 shows that 
the filamentary network identified lies only in the very 
densest part of the cloud, where the estimated mass per 
unit length value is signficantly above the thermal critical 
value (white contours in Figure I). We will return to this 
point further in Sections 3.2 and 3.3. 

2.3. Calculating Radial Column Density Profiles 

Once the filaments are identified, we measure the radial 
column density profiles along them. Since the filaments 
tend to converge toward the simulation centre, and some¬ 
times even intersect, care is needed to properly calculate 
the radial column density profile. First, we assign ev¬ 
ery pixel to the filament which it is closest to. Next, we 
exclude pixels which lie very close (< 0.01 pc) to two 
or more filaments - this value was chosen to provide a 
balance between not including too many locations which 
might provide non-representative measures of a given fil¬ 
ament profile, and not excluding too large a fraction of 
material around the filaments. We then calculate the 
mean column density of pixels in separation bins equal 
to the pixel size (^ 0.002 pc). Finally, to ensure that 
the filament profiles are accurate, we exclude the mea¬ 
surement for any radial bin where at least 25% of the 
total length of the filament, at that separation, was not 
included in the profile calculation. This final criterion 
ensures that all radial column density profile measure¬ 
ments used in our analysis are reliable - there are no 
cases where data from only a few pixels are used to in¬ 
fer the filament’s properties. We note that the above 
restrictions limit our analysis to a smaller range in radii 
than used in Smith et al. (2014), although the range is 
closer to Arzoumanian et al. (2011). Smith et al. (2014) 
analyze only the brightest one or two filaments in any 
given simulation snapshot, which ensures that the con¬ 
tamination in filament profiles will be minimal; with our 
inclusion of fainter filaments, only smaller radial separa¬ 
tions from a given filament spine are free from material 
from neighbouring filaments. Figures 3 and 4 show sev¬ 
eral example radial column density profiles which will be 
further discussed in Section 4. 

3. BASIC FILAMENT PROPERTIES 

Visual inspection of the resulting radial column den¬ 
sity profiles (e.g.. Figures 3 and 4) reveals a variety of 
characteristics. We expect that after the first turbulent 
shocks form a filamentary structure, gravity acts to con¬ 
tinue to concentrate mass onto these filaments, leading to 
higher and narrower peaks with time. An initially higher 
mean density should increase gravity’s pull and lead to 
a faster filament evolution. The presence of a magnetic 
field should cushion the initial turbulent compressions, 
reducing the amount of material initially in the filament, 
and giving the appearance of ‘fluffier’ filaments. The 
subsequent evolution of MHD filaments should then be 
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Fig. 1.— A comparison of the column density distribution projected along the X axis for the simulations. The top two rows show the 
500 M 0 simulations at 0.05, 0.1, and 0.15 Myr after the start of the simulation (left to right) for the MHD (top) and HD (middle) runs. 
The bottom row shows the MHD and HD 2000 Mq simulations at 0.05 Myr (left and middle) and the HD 500 Mq simulation at 0.2 Myr 
(right). All simulations are shown cropped to the inner 1.5 pc to better show the smaller-scale structure that forms. The greyscale range 
is the same in all panels, going from 0.01 to 10 g cm“^ (~ 3 x 10^^ to 3 x 10^^ cm“^) from black to white, with a logarithmic scaling 
applied. The overlaid contours show column densities of 0.02, 0.04, 0.075, and 0.2 g cm“^ (5.3, 1.1, 2.1, and 53 xlO^^ cm“^) in grey, white, 
blue, and red, respectively. If several assumptions are made, including that all pixels belong to cylindrical structures with a characteristic 
width of 0.1 pc, then the contours also correspond to mass per unit length values 0.5, 1, 2, and 5 times the critical mass per unit length 
at a temperature of 10 K (18 Mq pc“^). Note that these assumptions are poor for regions not associated with filamentary structure. See 
Section 3.4 for more detail. 
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Fig. 2. — Examples of the filamentary structure identified in the simulations. The top panel shows the filaments identified using DisPerSE 
in the X projection of the 500 Mq HD (left) and MHD (right) simulations at 0.15 Myr, while the bottom panel similarly shows the 2000 Mq 
HD (left) and MHD (right) simulations at 0.05 Myr. Each coloured line indicates a unique filament spine identified in that projection. 
Note that the top and bottom panels zoom in to different extents to best illustrate the central filamentary structure; similarly, each row 
has a different greyscale scaling applied - see the scalebar on the right hand side. In all panels, sink particles formed at the specified time 
are shown by the white stars; in all cases, their formation is confined to the central clustered part of the simulation. 
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Fig. 3.— Time evolution of the radial column density profile of a 
filament identified in the Z projection of the 500 Mq simulations. 
Solid lines show the reliable portion of the column density profile, 
while the dotted lines indicate less reliable measurements (i.e., data 
over less than 75% of the length of the filament was included). 
The error bars indicate the standard deviation of column density 
values at that radial separation. The MHD error bars are slightly 
shifted for better legibility. In this example, the filament continues 
to contract in both the MHD and HD simulations, although at a 
faster rate in the HD case. 

slowed relative to the HD case by gravity’s weaker pull 
on the the initial lower concentration of mass, and possi¬ 
bly also further action by the magnetic field, depending 
on its orientation. 

Broadly, these behaviours do hold - the visual impres¬ 
sion from watching movies of each simulation suggest 
this behaviour, nevertheless, we find instances of fila¬ 
ments dissipating over time, suggesting gravity was insuf¬ 
ficient to prevent the initial turbulent compression from 
re-expansion. In some of these cases, magnetic fields ap¬ 
pear to help to slow or prevent this re-expansion, causing 
the HD filament to have a higher and narrower peaked 



Fig. 4.— Time evolution of the radial column denstiy profile 
identified in the Z projection of the 500 MQsimulations. See Fig¬ 
ure 3 for details on the plotting conventions used. In this example, 
the filament contracts and then expands in both the HD and MHD 
simulations. 

profile than in the MHD case. In other instances, data 
excluded for one or more of the reasons mentioned above 
(difficulty in tracing the filament, or exclusion due to 
unreliability) also prevents the full influence of time or 
magnetic fields to be fully assessed. 

Despite this more complex behaviour, there are still 
several simple measures that we can make to gain insight 
into the behaviour observed. 

3.1. Filament Widths 

The conceptually simplest measureable filament prop¬ 
erty is its width. Although filament widths measured 
with Herschel span at least a factor of five (e.g.. Figure 7 
in Arzoumanian et al. 2011), it is often stated that fila¬ 
ments have a constant width of ^0.1 pc. Note, however, 
that Juvela et al. (2012a) find a larger scatter in filament 





























































FWHM values in their analysis of (different) Herschel 
data, although some of their filaments are much more 
massive and / or more distant than the Arzoumanian 
et ah (2011) sample. We measure the width of all of the 
filaments tracked in our simulations in the simplest pos¬ 
sible method - the extent of the radial profile at half of 
the peak value, i.e., the FWHM. Table 2 shows our re¬ 
sults, separated by time step, magnetic field, and mass. 
Included is the mean and standard deviation of FWHM 
values measured, along with the number of FWHM val¬ 
ues considered. Some filaments did not have reliable ra¬ 
dial column density profiles out to sufficiently large radial 
separations to allow the FWHM to be measured; these 
were excluded from the values given in Table 2. 

Although the dispersion is large, owing in part to the 
disparate behaviours discussed earlier, it is clear that on 
average, the filaments do behave as expected. Filaments 
generally get narrower with time and in higher mean den¬ 
sity environments, and are wider when magnetic fields 
are present. Furthermore, this trend is somewhat sub¬ 
tle: all of the simulation snapshots give filaments that 
have widths within the range that observers find. Heitsch 
(2013a,b) note that in accreting filament models, the fil¬ 
ament width can remain relatively constant throughout 
much of a filament’s evolution, if either the ram pressure 
from accreting material is small Heitsch (2013a) or in 
the case of a weakly magnetized accretion. Hennebelle & 
Andre (2013) propose a model wherein ion-neutral fric¬ 
tion dominating the dissipation of turbulence accounts 
for a relatively constant filament width of ^ 0.1 pc while 
Gomez & Vazquez-Semadeni (2014) suggest a constant 
width may be caused by a balance between large-scale ac¬ 
cretion onto filaments and accretion from the filaments 
onto the dense cores and stars forming within them. Our 
simulations do not contain ambipolar diffusive effects for 
the magnetic field, so that the substructure formed in¬ 
volves a balance between accretion, gravity, and turbu¬ 
lence. Future work will measure the accretion onto fila¬ 
ments versus the dense cores. 

We also tried fitting Gaussians to the filament profiles, 
with a constant background level set as a free parame¬ 
ter (fits not shown), similar to other work (Arzoumanian 
et al. 2011; Smith et al. 2014)^. Table 3 shows the equiv¬ 
alent FWHMs derived from the Gaussian fits for the fil¬ 
aments. Although the FWHM values tend to be smaller 
when using a Gaussian fit with a non-zero background, 
especially at earlier times (where the relative amplitude 
of the background is larger), the same general trends 
hold true: the filament width decreases with time, and 
tends to be larger when magnetic fields are present. The 
mean widths are still generally consistent with the ob¬ 
servations, although the range of widths is larger in the 
simulations, similar to the findings of Smith et al. (2014). 

3.1.1. Biases 

Smith et al. (2014) provide a detailed consideration of 
factors which can impact the measured filament width. 
For example, they find that when performing a Gaussian 
fit, the best-fit FWHM strongly depends on the radial ex¬ 


tent of the profile: including measurements from larger 
separations from the filament spine tends to increase the 
FWHM. Presumably this is at least partly caused by a 
lower background column density being fit for profiles 
that extend further from the filament spine. Our analy¬ 
sis tends to use a smaller radial extent than in Arzouma¬ 
nian et al. (2011) and especially Smith et al. (2014) due 
to potential contamination from other nearby filaments, 
which may explain why we tend to measure narrower 
filament widths in our Gaussian fits^. Measuring the 
FWHM directly from the profile will be robust to radial 
range variations, but has its own bias. Unresolved cen¬ 
tral filament peaks become lower with poorer resolution, 
which would change the peak column density used to es¬ 
timate the FWHM (R. Smith, priv. comm.). Both the 
biases in the FWHM and Gaussian-fitted width measure¬ 
ments are primarily systematic, affecting absolute rather 
than relative values. (We emphasize that this statement 
does not imply that the range in widths is invariant, but 
that the relative rankings likely are, i.e., the widest fila¬ 
ments appear to be the widest with any measure.) We 
expect then that our conclusions about the weak trends 
in width are therefore robust, a point supported by the 
similarity in behaviour using either width measurement. 

Finally, we note that the timescale over which we are 
able to analyze the filaments is relatively short: 0.2 to 0.3 
times the global free-fall time for the 500 Mq simulations, 
and 0.13 times the global free-fall time for the 2000 Mq 
simulations. Since the filaments form in the denser parts 
of the simulation, a larger number of local free-fall times 
would have elapsed. Nonetheless, analysis over a longer 
timescale could reveal stronger signs of filament evolution 
than we are able to probe here. 


3.2. Mass per Unit Length 

The simplest equilibrium model for a filament is that 
of the isothermal cylinder, presented in Ostriker (1964), 
where gravity is balanced by thermal pressure along an 
infinite cylinder. In this model, the stability of the cylin¬ 
der is controlled by the mass per unit length, Mune- The 
critical mass per unit length, in turn, depends only on 
the temperature: 


_ 2kBT ^ ^ 

lamnG G 


( 2 ) 


where Cg is the sound speed and G the gravitational con¬ 
stant (Ostriker 1964). Furthermore, Inutsuka & Miyama 
(1997) showed that isothermal filaments are unstable to 
axisymmetric perturbations of wavelength greater than 
about 2 times the filament diameter if the mass per unit 
length is close to this critical value. 

In our simulations, the temperature is set at a constant 
10 K, which implies = 18 Mq pc“^. Turbulent mo¬ 

tions can also provide additional support through raising 
the typical velocity dispersion of the gas above the ther¬ 
mal value; Heitsch (2013a) points out that non-thermal 
motions can be driven by the accretion of material onto 
the filament itself (see also Peretto et al. 2014). Fiege 
& Pudritz (2000) show that a more appropriate critical 


^ All of our fits (including those in Section 4) make use 
of the IDL mpfit routine by Markwardt (2009), available at 
http://www.physics.wisc.edu/ 

~craigm/idl/fitting.html 


® Smith et al. (2014) focus their analysis on the brightest one 
or two filaments in each of their simulations, which ensures that 
the contamination from other filamentary structures will be more 
minimal, even with a larger radial extent. 
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TABLE 2 

Filament FWHM values 


Mass 

(Mq) 

Time 

(Myr) 

HD - FWHM stats^ 
mean(pc) stddev(pc) 

MHD - FWHM stats^ 
mean(pc) stddev(pc) 

HD - FWHM stats'’ 
mean(pc) stddev(pc) 

MHD - FWHM stats'’ 
mean(pc) stddev(pc) 

500 

0.05 

0.212 

0.132 

0.262 

0.164 

0.211 

0.127 

0.268 

0.175 

500 

0.10 

0.105 

0.086 

0.176 

0.130 

0.066 

0.051 

0.116 

0.087 

500 

0.15 

0.079 

0.073 

0.130 

0.090 

0.050 

0.028 

0.104 

0.051 

500 

0.20 

0.058 

0.047 

N/A 

N/A 

0.044 

0.025 

N/A 

N/A 

2000 

0.05 

0.119 

0.113 

0.168 

0.147 

0.132 

0.137 

0.170 

0.142 


“ All measureable filaments included. ^ Only filaments where the FWHM could be measured at all times in HD and MHD were included. 

Note that the 500 Mq and 2000 Mq filaments are different samples. 


TABLE 3 

Filament Gaussian-fit FWHM values 


Mass 

(Mo) 

Time 

(Myr) 

HD - FWHM stats" 
mean(pc) stddev(pc) 

MHD - FWHM stats" 
mean(pc) stddev(pc) 

HD - FWHM stats'’ 
mean(pc) stddev(pc) 

MHD - FWHM stats'" 
mean(pc) stddev(pc) 

500 

0.05 

0.08 

0.07 

0.10 

0.07 

0.11 

0.10 

0.12 

0.09 

500 

0.10 

0.05 

0.03 

0.09 

0.05 

0.06 

0.03 

0.12 

0.06 

500 

0.15 

0.04 

0.02 

0.06 

0.03 

0.05 

0.02 

0.07 

0.03 

500 

0.20 

0.05 

0.02 

N/A 

N/A 

0.05 

0.02 

N/A 

N/A 

2000 

0.05 

0.04 

0.03 

0.10 

0.11 

0.05 

0.03 

0.15 

0.15 


All measureable filaments included. ^ Only filaments where a Gaussian could be fit at all times in HD and MHD were included. Note 

that the 500 Mq and 2000 Mq filaments are different samples. 


mass per unit length value is given by 

(3) 

where (a^) is the velocity dispersion including both ther¬ 
mal and non-thermal components. A careful analysis of 
the velocity fields would be required to determine pre¬ 
cisely how much nonthermal support is provided on the 
scales of interest; the approximation often assumed is 

(7^ = X (1 + A4 V3) (4) 

(e.g., Klessen et ah 2000). With A4 ^ 6 in our simu¬ 
lations, that would lead to raising by a factor of 

roughly 37, giving 670 Mq pc“^. 

In the case of a magnetized turbulent cloud, there is a 
magnetic correction that must be made to eqn 3. In the 
case that there is only a poloidal magnetic field and no 
toroidal (wrapped field) component, the magnetic pres¬ 
sure helps support the filament against gravity. Tak¬ 
ing eqn 27 of Fiege & Pudritz (2000) with F^ = 0 (no 
toroidal field), and using Fiege & Pudritz (2000) eqn 23 
to convert between the vertical magnetic field flux and 
magnetic field strength, the critical mass per unit length 
becomes 

= MZZ X (1 + Ml/2) (5) 

where is the turbulent Alfven Mach number 

whose Alfven speed, va, is B/^JAlTp. For super-Alfvenic 
turbulence, Ma < 1, this correction is slight. In our sim¬ 
ulation, however, Ma is 2.1 — 2.2 (see Table 1). Thus, 
our MHD turbulence is somewhat sub-Alfvenic, i.e., the 
magnetic field strength dominates the turbulence, and 
this implies that the critical line mass for our models is 
somewhat greater than the purely hydrodynamic turbu¬ 
lent case; ^ This result predicts that 

our MHD case should be considerably less susceptible to 
fragmentation than the hydro case. We note that as the 
turbulence is damped and the line width reduces to the 
thermal value, the relative magnetic contribution to the 


(thermal) line mass can become significantly more im¬ 
portant depending on the orientation of the field across 
the filament. 

There has been little time in these simulations for the 
turbulence to decay, but some of the power of the turbu¬ 
lence is on larger scales than our filaments, so the tur¬ 
bulent critical mass per unit length of 670 Mq pc“^ is 
an upper limit to the true effective critical mass per unit 
length of the simulated filaments. Indeed, the observed 
velocity dispersion in filamentary gas tends to be only 
of order twice the sound speed in nearby filaments that 
have peak column densities similar to those formed in 
our simulations (e.g., Hacar & Tafalla 2011; Hacar et al. 
2013; Arzoumanian et al. 2013; Kirk et al. 2013). Arzou¬ 
manian et al. (2013) furthermore find that for filaments 
with masses per unit length much higher than the criti¬ 
cal thermal mass per unit length, the velocity dispersion 
increases with increasing mass per unit length, a trait 
they attribute to infall of material onto the filaments. 
In our simulations, most if not all of of the turbulent 
motion is likely due to the remnants of the initial turbu¬ 
lence, given there has been little time for that to decay, or 
for infall to generate additional turbulent motions. As¬ 
suming the total filament velocity dispersion is twice the 
thermal value would increase by a factor of four. 

The estimated mass per unit length contours shown on 
Figure 1 illustrate that most of the dense filamentary 
structure is found in areas with above 4-5 for the 

500 Mq simulations, and even higher in the 2000 Mq 
simulations. Peretto et al. (2014) estimated that non¬ 
thermal support in SDC13 would contribute to lowering 
the mass per unit length values in the filaments to 1-2 
times the critical value from 4-8 times the critical value 
if non-thermal motions were not accounted for. Note 
that depending on the orientation, magnetic fields could 
either aid or hinder gravitational collapse. 

The degree of gravitational fragmentation of the fila¬ 
ments can, to some degree, be ascertained by the num¬ 
ber of sink particles that form in the simulations. It 
is notable that the sinks are typically found in or near 
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the filaments. It is also clear that the number of sinks 
that form in MHD simulations is smaller, sometimes no¬ 
tably so, in comparison with their hydrodynamic coun¬ 
terparts. Thus, using the data in Table 1, we see that in 
the 500 Mq simulation, while 16 sinks particles appear in 
the HD run, only 6 are apparent in the MHD case. The 
suppression is even greater in the 2000 Mq simulation 
(although in that case the runtime was shorter) where 
45 formed in the HD case as compared to 3 in the MHD 
case. Clearly, magnetic support is significantly reducing 
fragmentation. 

Non-isothermal equilibrium cylinder models have also 
been investigated (see Recchi et al. 2013, and the dis¬ 
cussion therein). In the case of a thermally-supported 
equilibrium cylinder, where the temperature gradually 
increases outward, similar to observations, Recchi et al. 
(2013) find that the mass per unit length which can be 
supported is only about 20 to 30% larger than in the 
isothermal case. Since the simulations we analyze are 
strictly isothermal, we do not consider this class of mod¬ 
els further here. 

3.3. Individual filament M/L measurements 

Measuring the total mass of each filament is difficult as 
can be seen in Figures 1 and 2, the filaments do not tend 
to have clearly defined outer boundaries. This challenge 
is exacerbated by the fact that many of the filaments 
he close to one another - the total mass cannot be de¬ 
rived by including material arbitrarily far away from the 
filament’s spine. We determined that the best way to 
estimate each filament’s mass was to include only mate¬ 
rial within the FWHM of the filament spine. While this 
will necessarily provide a lower limit to the true filament 
mass, lower and wider thresholds, such as the full width 
at quarter maximum, cannot be determined for too large 
a fraction of the total filament population (see discussion 
above). Using a constant width for the mass determi¬ 
nation would bias the estimates toward relatively lower 
values for the wider / fluffier filaments. We note that in¬ 
stead adopting a filament width based on the Gaussian 
fits discussed in Section 3 yields a similar behaviour to 
that discussed below. 

Using the estimated filament mass per unit lengths, we 
can test whether this metric is a useful predictor of fila¬ 
ment stability. A very simple proxy for the evolutionary 
path of a filament is to compare the peak column density 
(in the radial column density profile) at two neighbour¬ 
ing time steps. If the filament is contracting or accreting, 
the second peak should be higher than the first, while the 
reverse would be true for an expanding filament.^ We 
would therefore expect that higher mass per unit length 
values (above the critical value) would correspond to a 
ratio in peak column densities of greater than one (for 
the later time divided by the earlier time). 

Figure 5 shows the mass per unit length of each fila¬ 
ment compared to the ratio of the peak column density 
at the subsequent and current time steps. The thermal 

® We verified this assumption by comparing the ratio of filament 
FWHM values at subsequent times and found very good correla¬ 
tion: over 86% of filaments interpreted as contracting or expand¬ 
ing based on their peak column density ratio at neighbouring time 
steps show the same signature in their FWHM ratio. Allowing for 
slight measurement uncertainties (5% error in the ratios) gives an 
agreement between the two ratio measurements of just over 95%. 


critical mass per unit length for a temperature of 10 K 
is ^ ISMq pc“^, below nearly all of our measurements. 
Non-thermal motions likely contribute some amount of 
support (Section 3.2). In Figure 5 we show the critical 
mass per unit length assuming the total velocity disper¬ 
sion is twice the thermal value, which gives a value of 
^72 Mq pc“^. While this is a rough approximation, it 
appears to denote the level above which no filaments are 
found to be expanding. Regardless of precisely where 
the effective critical mass per unit length is drawn, we 
note a surprising result: many points occupy the bottom 
right quadrant, contracting filaments whose current mass 
per unit length implies gravitational stability. Similarly, 
there are several filaments whose mass per unit length 
ratio suggests graviational instability which are instead 
expanding. Although some of the presently contract¬ 
ing, low mass per unit length filaments (bottom right) 
may expand at time steps beyond what we can trace, 
the figure highlights the fact that predictions about the 
future evolution of filaments are incomplete when only 
column density information is available. This is also ap¬ 
parent in Figure 6, where the filament with the initially 
higher peak column density is the one which later re¬ 
expands. Although the individual mass per unit length 
values are not a good predictor of future evolution, the 
fact that there is a weak correlation between mass per 
unit length and peak column density ratio suggests that 
some (limited) insight into the bulk behaviour of fila¬ 
ments can be gained from the simple isothermal-with- 
turbulence model. 

3.4. Region-wide Mass per unit Length 

Finally, we can also get a rough idea of the stability of 
all the material in the simulation. The Hersehel Gould 
Belt team has provided an estimate of the mass per unit 
length at every pixel, for material in their curvelet map, 
i.e., that associated with structures having long axis ra¬ 
tios (e.g., Andre et al. 2010). These mass per unit length 
estimates are made under the assumption that each pixel 
is part of a filament or cylindrical structure with a typ¬ 
ical width of 0.1 pc. The contours in Figure 1 can be 
interpreted under a similar set of assumptions, although 
we emphasize that since our calculation includes the en¬ 
tire mass in the simulation, the equivalent mass per unit 
length values should not be over-interpreted in regions 
not associated with filamentary structure. The thermal 
critical value assumed for the contours in Figure 1 is that 
for a 10 K medium, i.e., 18 Mq pc“^. Most of 

the mass in filaments in the simulation lies above the 
critical mass per unit length value; if thermal pressure 
was the only source of support preventing gravitational 
collapse, we would expect to And sink particles forming 
throughout the simulation. Instead, all of the sink parti¬ 
cles form at mass per unit length values of at least 5 (red 
contour in Figure 1). Including nonthermal motions, as 
discussed earlier, raises the critical mass per unit length, 
therefore decreasing the ratio of the mass per unit length 
to the critical value. The mass per unit length of individ¬ 
ual filaments above which only contracting filaments are 
found, as discussed in Section 3.3, was a similar factor 
(4) above the thermal critical value. It therefore appears 
that only the densest parts of the simulation, where fil¬ 
aments are present, are likely to be sufficiently dense for 
gravitational collapse to occur. 
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Fig. 5. — The mass per unit length measured for a filament at a 
given time step versus the ratio in peak column densities for the 
subsequent versus current time step. Filaments which are contract¬ 
ing would be expected to have a ratio in peaks greater than one 
(right of the vertical dotted line), whereas filaments which are re¬ 
expanding would have a ratio in peaks of less than one (left of the 
vertical line). The thick horizontal dashed line indicates the esti¬ 
mated effective mass per unit length (assuming equal thermal and 
non-thermal contributions); filaments with values above this line 
would be expected to be contracting due to gravity, whereas those 
with lower values are stable against gravitational collapse. The thin 
horizontal dotted line indicates the critical mass per unit length 
assuming just thermal support; nearly all of the filaments have es¬ 
timated mass per unit length values well in excess of the thermal 
critical value. Note that the estimated mass per unit length values 
are all lower limits (see text for details). 

4. MODEL COMPARISONS 

We next compare the radial column density profiles 
obtained to several cylindrical equilibrium models: the 
isothermal cylinder, modified isothermal cylinder, and 
pressure confined isothermal cylinder, described in de¬ 
tail below. In the simplest model, the isothermal cylin¬ 
der, thermal pressure balances gravity along an infinite 
cylinder, leading to a 3D density profile decreasing as 
at large radii, or the column density varying as r~^ 
(Ostriker 1964). Herschel teams have found that the 
observed column density profiles are shallower than the 
isothermal cylinder model and and instead use a mod¬ 
ified profile, also referred to as a ‘Plummer-like’ profile 
(Nutter et al. 2008; Smith et al. 2014; Plummer 1911), 
where the power law exponent is an additional fitted pa¬ 
rameter: 

S(r)=^p- ^ (6) 

{l + {r/Rfir) - 

Here, H is the column density, pc is the central density, 
Rfi represents the scale of the inner flat portion of the 
profile, p is the power law index (with a value of 4 for the 
original Ostriker model), and Ap is a geometrical factor 
given by: 

A = ^ f ( 7 ) 

^ cosi (1 + 


where i is the (unknown) inclination of the filament on 
the plane of the sky, assumed to be 0 (Arzoumanian et al. 
2011). The best-fitting value of p often tends to range 
between 1.5 and 2.5, a much shallower drop-off than in 
the p = 4 Ostriker (1964) model (e.g., Alves et al. 1998; 
Lada et al. 1999; Arzoumanian et al. 2011; Malinen et al. 
2012; Juvela et al. 2012b; Hill et al. 2012; Palmeirim et al. 
2013). Some filaments, however, have been observed with 
column density profiles which are consistent with a p = 4 
isothermal model, (e.g. Nutter et al. 2008; Pineda et al. 
2011; Hacar & Tafalla 2011; Bourke et al. 2012), while 
Contreras et al. (2013) find p = 4 provides a good fit 
around star-forming clumps and p = 2 is better in the 
inter-clump areas. Theoretically, shallow radial column 
density profiles are consistent with equilibrium isother¬ 
mal cylinder models that include helical magnetic fields 
(Fiege & Pudritz 2000). Smith et al. (2014) show that 
p ^ 2 profiles are the norm for prominent filaments in 
hydrodynamic simulations without magnetic fields, re¬ 
gardless of the type of turbulence considered. 

We also applied the equilibrium model of Fischera & 
Martin (2012), in which pressure from the medium sur¬ 
rounding the filament is also included in the force bal¬ 
ance. In this analytic formulation, the two quantities of 
interest are P, the pressure from the external medium, 
and /, the ratio of the mass per unit length to the crit¬ 
ically stable value for the Ostriker cylinder. The full 
profile is given by: 

i-/~A 

(^v'/(i-/)(i-x’)+^/FkA) 

where Nh is the column density in number units, x is a 
scaled radial coordinate, G is the gravitational constant, 
p is the mean molecular weight, and tuh is the mass of 
a hydrogen atom (Fischera & Martin 2012). The main 
effect of the pressure, P is on the height of the central 
column density peak, while / controls the shape of the 
profile (Fischera & Martin 2012). [The temperature of 
the gas is also fit as part of the scale factor converting 
the radial coordinate x into a physical radial separation.] 
In their re-analysis of the Herschel filaments in Polaris, 
IC 5146, and Aquila, Fischera & Martin (2012) find that 
their pressure equilibrium model also provides a good fit 
to the filaments, with external pressures consistent with 
the range expected in the ISM. 

In all of the models, a non-zero background column 
density can be included as a free parameter, i.e., 

E(r) = Y,modei{r) + Eq (9) 

where Sq is a constant. 

Most if not all of the filament analyses include a back¬ 
ground column density term (e.g., Arzoumanian et al. 
2011; Juvela et al. 2012b; Fischera & Martin 2012; 
Palmeirim et al. 2013; Smith et al. 2014); Herschel anal¬ 
yses in fact allow the background to be fit by a linearly 
varying background column density: S(r) = T^modei{p) + 
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i;o + Hi(r) (see Appendix B of Palmeirim et al. 2013), al¬ 
though the fits are generally similar when just a constant 
background column density is adopted (D. Arzoumanian, 
priv. comm.). We tried fitting the profiles with and with¬ 
out a background column density and found that includ¬ 
ing the background generally produced superior fits. In 
the case of the modified isothermal cylinder model (eqn 
5), including a background term decreased the central 
density, and allowed for a narrower peak to be fit, better 
matching the filament profiles. The difference was most 
pronounced for the case of the pure isothermal cylinder 
model (eqn 5 with p = 4), where few cases produced good 
fits without a background column density. The pressure 
confined model (eqn 7) almost never converged to a sat¬ 
isfactory fit without a background column density term. 
Appendix A shows the results of fits with no background 
column density included. 

Figure 6 shows the best fit models for one example 
filament at 0.1 Myr. Following the general behaviour 
seen in the simulations, the filament in the MHD simu¬ 
lation is ‘fluffier’ (wider and lower peak column density) 
than in the HD simulation. This is a consequence of the 
significant amount of magnetic pressure support of the 
filaments as noted earlier. 

4.1. Modified Isothermal Cylinder Model Fits 

As can be seen in Figure 6, the isothermal and modi¬ 
fied isothermal cylinder model provide near-identical fits, 
when a background column density term is included. In 
both the HD and MHD profiles shown, the models have 
peak values within the errors of the simulated radial pro¬ 
file. As shown in Appendix A, the models differ by a 
greater amount when the background column density is 
fixed to zero, with the pure isothermal model then gen¬ 
erally providing a very poor fit. 

Table 4 gives the median model fit parameters for fila¬ 
ments where a fit was possible at every time step in both 
the HD and MHD simulations, to allow for any trends 
in the time evolution to be followed. The typical power- 
law slopes found (p ^ 1.3 — 2) are similar to the best-fit 
values in Arzoumanian et al. (2011). The typical disper¬ 
sion (standard deviation) in the model fit parameters is 
of the same size as the median values given in the ta¬ 
ble; in all cases, there is no clear trend in the best fit 
parameters evolving as a function of time. Juvela et al. 
(2012a) and Smith et al. (2014) point out that the modi¬ 
fied isothermal fit is partially degenerate between the fit 
parameters, which could hide evolutionary trends. Smith 
et al. (2014) demonstrated that while no time evolution 
was seen in their best-fit parameters when all were free 
to vary, fixing Rfi gave best-fit central densities which 
clearly increased with time. We performed the same test 
and found a similar result, as shown in Table 5, although 
the combination of expanding and contracting filaments 
at later time steps diminishes the strength of the evolu¬ 
tionary trends. 

Although the variation between filaments in a single 
simulation snapshot is larger than any general change 
in filament behaviour as a function of time or initial 
conditions, the model fits are consistent with the gen¬ 
eral behaviours noted earlier. Magnetic fields produce 
puffier (lower pc and higher Rfi) filaments, and an ini¬ 
tially higher density tends to lead to more peaky fila¬ 
ments (higher pc and lower Rfi)] all of these differences 
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Fig. 6. — The radial column density profile for the filament in 
Figure 3 at 0.1 Myr, including the best-fit models with a back¬ 
ground column density term. The top panel shows the profile and 
fits for the HD simulation, as well as the residuals between the 
models and simulation, while the bottom panel shows the same for 
the MHD simulation. The solid black line indicates the mean col¬ 
umn density at each radial separation, while the error bars denote 
the standard deviation in values at various radial separations. The 
model lines shown are the purely isothermal cylinder (darkest, dot¬ 
ted line), the modified isothermal cylinder (less dark, dashed line), 
and the pressure confined cylinder (lightest, dash-dotted line). 

are smaller than the scatter in best fit values for a given 
simulation snapshot, and would require tracking over a 
longer time period to better measure time evolution. 

As Arzoumanian et al. (2011) found, we also find that 
the isothermal cylinder model tends to provide a worse 
fit to the radial column density profile than the modified 
isothermal model, although this is much more noticeable 
in the case of a zero background column density (see Ap¬ 
pendix A). Typically, the best fit power law index, p is 
between 1.3 and 2.0, within the range found by Arzouma¬ 
nian et al. (2011) and others, and much less than p = 4 
for the pure isothermal model. We also find pc is around 
10^ to 10^ cm“^, and Rfi ranges between roughly 0.001 
to 0.1 pc, consistent with Juvela et al. (2012b). The re¬ 
lationship between Rfi and the FWHM is dependent on 
the power law of the profile, and that a smaller value of 
Rfi is expected based on the values of p fitted. Explic¬ 
itly, FWHM = 2 X — IRfi, or 3A6Rfi when 

p = 2. 

4.2. Pressure-Confined Isothermal Cylinder Model Fits 
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TABLE 4 

Best fit parameters for modified isothermal cylinder, with 

BACKGROUND 

Mass Time HD « MHD « 


(Mq) 

(Myr) 

pc 

Rfi 

P 

So 

pc 

Rfi 

P 

So 

500 

0.05 

0.9 

3.0 

2.3 

5.4 

1.4 

3.0 

2.4 

7.7 

500 

0.10 

14.3 

0.6 

1.5 

3.2 

5.5 

2.0 

1.8 

3.5 

500 

0.15 

8.8 

0.8 

1.8 

8.7 

1.7 

3.6 

3.1 

4.5 

500 

0.20 

4.5 

2.3 

4.2 

9.7 

N/A 

N/A 

N/A 

N/A 

2000 

0.05 

8.4 

1.0 

1.8 

29 

6.7 

1.1 

1.4 

23 

500 

alF 

4.6 

1.7 

2.0 

4.6 

1.7 

2.7 

2.0 

4.3 

2000 

alF 

8.4 

1.0 

1.8 

29 

6.7 

1.2 

1.4 

23 


“ Mean of the best fit values for filaments fit at all times, with a 
non-zero background allowed for the fit: the central density, pc 
(in units of 1(D cm“^), the central flat radius, Rfi (in units of 
0.01 pc), the exponent, p, and Sq, the background column 
density term (in units of 10^^ cm“^); the standard deviation is 
often comparable in magnitude to the mean, with values of 
0.6-18, 0.4-2.4, and 0.2-3.5 for pc, Rfi and p respectively, in the 
same units. The background column density tends to have 
standard deviations larger than the mean value, and the mean is 
usually significantly larger than the median. 

^ Values of all profiles where a fit was possible are included here, 
i.e., relaxing the requirement of a fit for both HD and MHD, and 
for the 500 Mq simulations, a fit at all time steps. 


TABLE 5 

Best fit parameters for modified isothermal cylinder, with 
BACKGROUND, Rfi FIXED AT 0.01 PC. 


Mass 

(Mo) 

Time 

(Myr) 

Pc 

HD « 

P 

So 

Pc 

MHD “ 

P 

So 

500 

0.05 

1.6 

1.6 

3.2 

1.6 

1.4 

2.5 

500 

0.10 

2.7 

1.8 

7.0 

1.9 

1.5 

2.6 

500 

0.15 

3.3 

2.1 

8.6 

2.2 

1.6 

4.2 

500 

0.20 

3.6 

2.0 

7.3 

N/A 

N/A 

N/A 

2000 

0.05 

6.4 

1.9 

56 

6.1 

2.5 

55 

500 

alF 

1.6 

1.6 

6.7 

1.6 

1.5 

2.8 

2000 

alF 

7.3 

2.3 

68 

8.8 

2.1 

58 


“ Mean of the best fit values for filaments fit at all times, with a 
non-zero background allowed for the fit and the central flat 
radius, Rfi fixed at 0.01 pc: the central density, pc (in units of 
10^ cm“^), the exponent, p, and Sq, the background column 
density term (in units of 10^^ cm“^); the standard deviation is 
usually slightly less than the mean, with values of 0.9-2.3, and 
0.2-0.7 for pc and p respectively, in the same units. The 
background column density tends to have higher standard 
deviations than the mean values, and the mean is often much 
higher than the median value. 

^ Values of all profiles where a fit was possible are included here, 
i.e., relaxing the requirement of a fit for both HD and MHD, and 
for the 500 Mq simulations, a fit at all time steps. 

The pressure confined isothermal cylinder model also 
generally provides a good fit to the filament profiles, al¬ 
though the HD example shown in Figure 6 does a poor 
job of capturing the peak column density. In most cases, 
the best-fit model is very similar to the best-fit modified 
isothermal cylinder model. The mean and standard de¬ 
viation temperature of all fits are 15 ± 14 K, with a small 
tail in the temperature distribution extending out to 
100 K; only 20% of the temperatures fit were above 20 K; 
the median temperature fit is 11 K. The typical external 
pressure fit was Pextj^B = 4 ± 3 X 10^ cm^ K ^ , with 
the tail in the distribution extending up to 10^ cm^ K“^. 
The typical shape parameters fit were / = 0.76 ±0.18 
(keeping in mind / must be between 0 and 1). The fitted 
temperature and external pressure values are physically 
reasonable - the temperatures are generally similar to 
the simulation’s constant 10 K, and the typical external 


pressure is is the same range as those fitted and esti¬ 
mated to be reasonable in nearby molecular clouds such 
as Perseus (Kirk et al. 2006). 

We searched for evolutionary trends within the fitted 
model parameters, but did not find any over the time pe¬ 
riod analyzed. The only discernable trend was that the 
external pressures fit tended to be higher in the 2000 Mq 
simulations than the 500 Mq simulations, with typical 
values of 6 ± 6^ cm^ and 1 ± 2 x 10^ cm^ re¬ 
spectively. Since the initial density of the 2000 Mq sim¬ 
ulation was higher, the external pressure caused by the 
weight of overlying material within the region would be 
expected to be higher. 

Finally, we made a general comparison of the goodness- 
of-fit of the various models, by comparing the typical 
(mean and standard deviation) chi-squared values of all 
fits. The standard deviation in chi-squared values is sev¬ 
eral times larger than the mean for any of the models 
fit, making a distinction in the overall goodness-of-fit 
between the models tenuous. Nevertheless, including 
a background column density term for the isothermal 
cylinder model made a substantial difference: the mean 
chi-squared value for fits with a background included is 
more than 3.5 times smaller than when the background 
is zero. The difference is much less pronounced for the 
modified isothermal model where including a background 
column density decreases the mean chi-squared value by 
40%. Allowing the power law to vary (purely isothermal 
model versus modified isothermal model) yields a 30% 
improvement in the mean chi-squared value, while the 
pressure confined model has a mean chi-squared value 
8% lower than the modified isothermal model. 

5. EFFECT OF RESOLUTION 

Finally, we examine the effect of resolution on our re¬ 
sults. Real observations of filamentary structures are 
complicated by both instrumental effects (including res¬ 
olution and system noise) and physical effects (how the 
flux emitted at a specific wavelength relates to the in¬ 
trinsic column density of material), which make direct 
comparisons between simulations and observations diffi¬ 
cult (e.g., Goodman et al. 2009). Here, we consider only 
the simplest factor, spatial resolution, to represent what 
‘perfect’ observations would be able to reveal. For the 
analysis presented here, we assume that the observed sys¬ 
tem is located 140 pc away, representing the very nearest 
molecular clouds and therefore also a best-case scenario. 
While we tested a variety of resolutions, we show results 
from three cases which are representative of the present 
single-dish facilities able to map large areas of the sky in a 
reasonable amount of time. The first case we consider is a 
resolution of 8". This corresponds to the resolution of the 
JCMT at 450 /im (Holland et al. 2013), and is also similar 
to the 1" resolution of the recent CLASSY survey (e.g., 
Fernandez-Lopez et al. 2014) which studies in detail three 
nearby molecular clouds in unprecedented detail with the 
CARMA interferometer. The second case we consider is 
a resolution of 18.2" corresponding to the resolution of 
Herschel column density maps using the group’s latest 
standard method (e.g., Palmeirim et al. 2013), with the 
resolution corresponding to that of their 250 /rm observa¬ 
tions. The third and final case we consider is a resolution 
of 36.9" corresponding to the resolution of Herschel at 
500 yum, and earlier Herschel column density maps (e.g. 
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Konyves et al. 2010; Arzoumanian et al. 2011). 

Figure 7 shows an example of the effect of a resolution 
on the column density in the simulation. As is clearly 
illustrated by this comparison, at the longest wavelength 
Herschel resolutions, much of the fine filamentary struc¬ 
ture is lost, while SCUBA2 at 450 yum / CLASSY does 
a better job. 

We also examined the quantitative effect of the resolu¬ 
tion on our results. To do this, we applied the same fila¬ 
ment definitions and re-ran our analysis. We attempted 
to correct for the resolution in a similar manner to real 
observational analyses: for the direct filament FWHM 
measurements, we deconvolved the values with the res¬ 
olution. For the column density profile models, we con¬ 
volved the model with the ‘beam profile’ / resolution be¬ 
fore fitting. We find that despite attempting to correct 
for the resolution during analysis using the standard ob¬ 
servational techniques, the poorest resolution still gives 
biased results. A full comparison of the effect of reso¬ 
lution on all of our measured quantities is given in Ta¬ 
ble 6: we calculated the mean and standard deviation 
of the ratio in values between the decreased resolution 
and original results. In summary, while there is a sig¬ 
nificant amount of scatter, typically the resolution only 
produces a somewhat noticeable effect for the 36.9" or 
Herschel 500 yum case. The effects tend to be in the direc¬ 
tion expected: poorer resolution leads to higher widths 
and lower central densities. We emphasize that these 
tests correspond to the best possible case for these in¬ 
struments. Very few star-forming clouds are as close 
as 140 pc; even Perseus is nearly twice as far away at 
^250 pc, and many other Gould Belt Survey clouds are 
at a similar or larger distance. Filaments surveyed in 
the Herschel Hi-Gal Survey are several times more dis¬ 
tant still. We ran additional tests (not shown) with even 
poorer spatial resolution, corresponding to filaments at 
these further distances, and found that the biases in ob¬ 
served quantities becomes much more noticeable in those 
cases. 

Smith et al. (2014) tested the effect of resolution on 
measured filament widths and find relatively little effect. 
Their test used a degraded resolution of 0.0086 pc (R. 
Smith, priv. comm.), corresponding to a resolution of 
^ 13" at 140 pc. Juvela et al. (2012b) also examined 
the effect of resolution on simulated filaments, using an 
angular resolution of 40" at 93, 186, and 371 pc. They 
found that at the larger distances, the filament central 
densities were the most biased, while the width, mass per 
unit length, and power law slope changed by less. Both 
of these results are consistent with our findings at similar 
resolutions. 

Beyond the implications to the measured filament 
properties presented here is the presence and charac¬ 
terization of substructure within the filaments. Aver¬ 
aging a radial column density profile along a filament 
hides much of the information on smaller-scale structure 
within the filaments in the simulations - see, for exam¬ 
ple, the comparison of 2D and 3D column density ver¬ 
sus density profiles of simulated filaments in Gomez & 
Vazquez-Semadeni (2014) and Smith et al. (2014). Both 
observations (Hacar et al. 2013; Henshaw et al. 2014; 
Fernandez-Lopez et al. 2014), and simulations (Moeckel 
& Burkert 2014; Smith et al. 2014) suggest that filaments 
may in fact be composed of multiple strands of dense 


TABLE 6 

Effect of Resolution on Eit Parameters" 


Resolution 

FWHM^ 



JCMT-450 lira 

1.2 ±0.4 

1.1 ±0.2 

1.1 ±0.2 

Herschel-250 [im 

1.3 ±0.5 

1.1 ±0.2 

1.2 ±0.3 

Herschel-500 iim 

1.6 ±1.2 

1.3 ±0.4 

1.4 ±0.6 

Resolution 

Pc 

Rh 

p" 

JCMT-450 lira 

0.8 ±0.3 

1.8 ±1.3 

1.2 ±0.4 

Herschel-250 /xm 

1.1 ± 1.1 

1.5 ±1.1 

1.1 ±0.5 

Herschel-500 fim 

0.5 ±0.3 

3.4 ±2.6 

1.6 ±0.7 

Resolution 

pi 

TDd 

jF 

JCMT-450 /im 

1.1 ±0.6 

1.1 ±0.3 

- 

Herschel-250 /xm 

1.1 ±0.6 

1.1 ±0.3 

- 

Herschel-500 /im 

0.8 ±0.4 

1.5 ±0.7 

- 

Resolution 


Pext 

fe 

cyl 

JCMT-450 lira 

1.1 ±0.2 

1.2 ±0.5 

1.0 ±0.2 

Herschel-250 /xm 

1.2 ±0.6 

1.3 ±0.6 

1.0 ±0.2 

Herschel-500 /rm 

1.4 ±0.7 

1.4 ±0.8 

1.0 ±0.3 


“ Mean and standard deviation in the ratio between decreased 
resolution fits and original values (all quantities). 

^ Filament FWHM widths, Gaussian-fitted widths, and mass per 
unit length ratio. 

Parameters from the modified isothermal cylinder model. 

^ Parameters from the isothermal cylinder model (power law 
exponent fixed). 

® Parameters from the pressure-confined isothermal cylinder 
model. 

gas, perhaps woven together, which are often difficult to 
identify separately with the present observational capa¬ 
bilities; higher spatial resolution and / or inclusion of the 
line of sight velocity are essential. 

In our hydrodynamic simulations too (see Figure 7), we 
see evidence of more complex structure on smaller scales, 
although a full 3D consideration of this is beyond the 
scope of the present paper. Nonetheless, our results cou¬ 
pled with new results emerging from observations such 
as Fernandez-Lopez et al. (2014) suggest that the key to 
a deeper understanding of filamentary structure requires 
high (spatial and velocity) resolution as well as more so¬ 
phisticated tools by which to model the structure. 

6. DISCUSSION & CONCLUSIONS 

We simulate the formation of filaments within a clump- 
scale volume, investigating the role of magnetic fields on 
the evolution of filaments column density. Starting with 
500 Mq (or 2000 Mq) of gas within a 2 pc cube, we track 
the evolution of filamentary structures over 0.15 Myr 
with and without the presence of a magnetic field. Tur¬ 
bulence remains strong throughout the simulations as 
there has been been insufficient time for it to damp sig¬ 
nificantly. 

These analyses provide an important complement to 
the recent filamentary analysis in simulations by Smith 
et al. (2014). Those authors investigated the effect of 
different initial modes of turbulence (e.g., solenoidal vs 
compressive), whereas our analysis examines the effect 
of magnetic fields and gravity (through varying the ini¬ 
tial mean density). Other more subtle differences are 
also important to note as well. Smith et al. (2014) as¬ 
sumed a uniform density initial sphere, surrounded by 
a warm diffuse medium, whereas we assume a sphere 
with a radially-decreasing density surrounded by a vac¬ 
uum. Girichidis et al. (2011) demonstrate that the ini¬ 
tial density distribution can have a marked effect on the 
large-scale structure which forms later in the simulation. 
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Fig. 7.— A comparison of the column density in the original simulation and ideal observations at 8.0'', 18.2", and 36.9" at 140 pc, 
corresponding to SCUBA2 at 450 fim (or CLASSY) and Herschel at 250 /rm and 500 fim respectively (left to right). The top panel shows 
the full 2 pc simulation box for the pure HD simulation with 500 Mq viewed along the x axis at 0.15 Myr. The bottom panel shows a 
zoomed in view to the box indicated in the top panel. The blue circles indicate the beamsize / resolution for each panel. 


Our simulations do not include the effect of radiation or 
simple chemistry, while Smith et al. (2014) does include 
both; we expect these effects to become more important 
at later times (once massive YSOs begin ionizing their 
natal environments) and when making synthetic obser¬ 
vations of molecular line emission (where the presence or 
absence of various molecular species has a large effect). 
The base numerical codes used are also different - Smith 
et al. (2014) use AREPO, a hybrid code, while we use 
FLASH , which is an adaptive mesh refinement based code. 
Despite these significant differences, both Smith et al. 
(2014) and our analyses do identify filaments which have 
properties broadly similar to observed filaments, which 
appears to speak to the universality of filament formation 
under a variety of conditions. 

We also note one major difference from Smith et al. 
(2014) in the analysis stage: Smith et al. (2014) focus 
their analysis on the brightest one or two filaments in 
each simulation whereas we also include fainter / less 
dense filaments in our analysis. In this respect, our anal¬ 
ysis gives a better direct comparison with observational 
surveys, where many filaments are identified in any given 
star-forming region. Our approach limits our analysis 
of the filament column density profiles to a smaller ra¬ 
dial separation from the filament spine, since the fainter 
filaments are more liable to have their profiles contami¬ 
nated by nearby non-related substructure than brighter 
filaments are (the filling factor is much larger for all faint 
plus bright filaments than it is for just bright filaments). 
On the other hand, inclusion of fainter filaments allows 
us greater sensitivity to less gravitationally bound fila¬ 
ments, which may be more liable to re-expand with time; 
such behaviour was not noted in Smith et al. (2014), pre¬ 
sumably because the dominant filament in each simula¬ 


tion will continue to contract and accrete new material 
throughout time. 

Our main findings are as follows: 

1. Magnetic fields have a strong influence on filamen¬ 
tary structure. Even with a mass to magnetic flux 
ratio which is supercritical by a factor of ^2, there 
are notable differences from the purely hydrody¬ 
namic simulation. Filaments formed in the mag¬ 
netic case tend to be wider, less centrally peaked, 
and evolve more slowly than filaments seen in the 
purely hydrodynamic case. These differences are 
most apparent through a visual comparison (e.g.. 
Figures 2, 3, and 4), and are less discernable in 
quantitative measures due to the large variation in 
filament properties at any given snapshot. 

2. The magnetic field can have a strong effect on the 
fragmentation of filaments. In our simulations, 
magnetic fields are able to significantly suppress 
the formation of cores, since its energy density ex¬ 
ceeds that of the turbulence. The turbulence is 
sub-Alfvenic {Ma = UA/cr^2.1 — 2.2) and so the 
magnetic field increases the critical turbulent line 
mass by a factor of 3.2 to 3.3. This accounts for 
the less condensed structure of the magnetized fil¬ 
aments, and their notably less fragmentation. 

3. The simulated filaments have properties consis¬ 
tent with observations. The radial column den¬ 
sity profiles of the filaments are well-described by a 
Plummer-like or modified isothermal cylinder pro¬ 
file. The power-law slope tends to be around 1.3 
- 2, similar to the range of 1.5 - 2.5 found in Her¬ 
schel data by Arzoumanian et al. (2011). The cen- 
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tral density tends to be of order 10^ cni“^ for the 
500 Mq simulations and closer to 10^ cm“^ for the 
2000 Mq simulations; the inner flat radius is a few 
hundredths of a parsec in both cases. The 500 Mq 
typical central densities, as well as the inner flat 
radii and power law slopes are consistent with those 
in Juvela et ah (2012a), also based on Herschel ob¬ 
servations. The pressure confined cylinder model 
of Fischera & Martin (2012) also provides a rea¬ 
sonable fit to the radial column density profiles, 
with typical temperatures, external pressures, and 
shape parameters fit of 15 K, 4 x 10^ cm^ K and 
0.76 for T, PextjkB^ and / respectively. 

4. Filaments have diverse evolutionary paths. At any 
given snapshot in time, the simulation reveals a va¬ 
riety of filaments. Some continue to radially con¬ 
tract and accrete material throughout the simu¬ 
lation, and will presumably continue on to form 
stars along their length. Other filaments halt 
in their contraction and expand into the ambient 
medium before the end of the simulation. Given 
the relatively short time duration of our simula¬ 
tions, we expect that even more diverse evolution¬ 
ary paths could be possible throughout the lifetime 
of a molecular cloud. 

5. The mass per unit length of a filament in a given 
snapshot provides only a weak diseriminant between 
the eontraeting and expanding filaments. Over 
most of the range of mass per unit length values, 
filaments can be either contracting or expanding. 
Above roughly 72 Msol pc“^, the critical value for 
filaments supported equally by thermal and non- 
thermal support, nearly all filaments are contract¬ 
ing. Velocity and magnetic field information are 
clearly required to determine the evolutionary state 
of a filament unambiguously. 

6. Turbulenee plays an important role in the mass per 
unit length of filaments. The filaments in which 
stars appear in our simulations have critical mass 
per unit lengths that are dominated by turbulent 
velocity dispersion and not just thermal values. 
This arises in these simulations since they are much 
less than a free-fall time old, so that turbulence 
has not had the opportunity to damp significantly. 
Hersehel results indicate that the thermal criti¬ 
cal value of Mifff is a good diagnostic of star¬ 


forming ability therefore suggests that turbulence 
is much weaker within these filaments, perhaps be¬ 
cause their natal clouds are more evolved. 

7. Filament widths are mildly influeneed by environ¬ 
ment. Filaments tend to be narrower at later times, 
when formed in higher density environments, or 
without the presence of magnetic fields, but these 
trends are all weak, given the mixture of contract¬ 
ing and expanding filaments at every snapshot in 
the simulations. Stronger trends in the evolution 
in filament properties as a function of time might 
become apparent with a longer timescale for analy¬ 
sis, especially if the filaments which re-expand be¬ 
come sufficiently diffuse to become undetectable at 
later times. We find the mean fiiament FWHM 
ranges from 0.06 to 0.26 pc across our simulations 
at varying times, with most being around 0.1 pc to 
0.15 pc^^. The range in filament FWHM values for 
any given simulation snapshot is, however, larger 
than the range seen in the analyses of Arzouma¬ 
nian et al. (2011). For filaments which are contract¬ 
ing, a combination of the decaying turbulence and 
gravity is likely responsible for the evolution. This 
may in part explain the relative constant filament 
widths discussed in Arzoumanian et al. (2011), al¬ 
though resolution may also be an important factor 
(Fernandez-Lopez et al. 2014), and observational 
biases and measurement methods could be playing 
a role (Heitsch 2013a; Smith et al. 2014). 

8. Filaments have eomplex struetures. The radial 
column density profiles of the filaments reveal a 
wealth of sub-structure, particularly in the pure 
hydrodynamic simulation where features tend to 
be sharper. Some of these substructures appear 
suggestively like the intertwined filament bundles 
found by Hacar et al. (2013) (see for example the 
bottom left panel in Figure 7); high-resolution ob¬ 
servations of a suite of filaments will be necessary 
to show how commonplace this phenomenon is. 
Other simulations, including Smith et al. (2014) 
and Moeckel & Burkert (2014) also find compiex 
3D fiiamentary substructure. 

Future work wiii inciude a 3D anaiysis of the fiia¬ 
ment s formed in these simulations, including their ve¬ 
locity structure and accretion rates, and the relationship 
with the magnetic field geometry. 


See Section 3.1 for a discussion on the biases and uncertainties 
in absolute filament widths in our analyses. 


APPENDIX 

Here, we examine the results of fits to the radiai coiumn density profiies of the fiiaments when the background 
coiumn density is forced to be zero (see discussion in Section 4). Figure 8 shows the best fit modeis for the isothermai 
and modified isothermai modei, for a fiiament identified in both the HD (top panei) and MHD (bottom panei) 500 Mq 
simuiations, to be contrasted with Figure 6 for the same fits where background coiumn density is aiiowed to be non¬ 
zero. Note that the best-fit pressure confined modei piotted in Figure 8 does inciude a background coiumn density 
term, and is identicai to the fit shown in Figure 6. Contrary to the case (Section 4) when the background coiumn 
density is fixed to zero, the three models differ from each other by a greater amount, and the pure isothermal model 
is generally a poor fit to the profiie. 

Tabie 7 gives the median modei fit parameters for the isothermal and modified isothermai profiies for aii fiiaments 
fitted at every time step in both the HD and MHD simuiations. Comparison of these fit vaiues with those given in 
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Fig. 8. — The radial column density profile and best fit models for a filament, including the fitting residuals where the background column 
density level for the isothermal and modified isothermal cylinder models is fixed to zero. See Figure 6 for the plotting conventions used. 
Note that the vertical range in the residual plots differs between the upper and lower profiles. In this case, the pure isothermal cylinder 
model is not a good match to the profile. 


TABLE 7 

Best fit parameters for modified isothermal cylinder, no background 


Mass 

(Mq) 

Time 

(Myr) 

Pc 

HD “ 

Rfi 

p 

Pc 

MHD “ 

Rfi 

P 

500 

0.05 

1.1 

1.7 

1.4 

0.7 

2.0 

1.4 

500 

0.10 

3.1 

0.7 

1.4 

0.7 

1.8 

1.4 

500 

0.15 

3.9 

0.6 

1.5 

2.3 

1.2 

1.4 

500 

0.20 

2.6 

0.8 

1.5 

N/A 

N/A 

N/A 

2000 

0.05 

11.1 

0.7 

1.3 

7.9 

0.6 

1.2 

500 

alP 

2.4 

0.7 

1.4 

1.0 

1.6 

1.4 

2000 

alP 

11.1 

0.5 

1.2 

10.7 

0.6 

1.3 


^ Median of the best fit values for filaments fit at all times: the central density, pc (in units of 10^ cm ^), the central flat radius, Rfi (in 
units of 0.01 pc), and the exponent, p; the standard deviation is often comparable in magnitude to the median, with values of 0.7-32, 

0.7-4.4, and 0.2-0.5 respectively, in the same units. 

^ Values of all profiles where a fit was possible are included here, i.e., relaxing the requirement of a fit for both HD and MHD, and for the 

500 Mq simulations, a fit at all time steps. 


Table 4 shows that the best fit profiles tend to have narrower peaks (better matching the filament profiles) when a 
background column density included, most of this increased narrowness is accounted for by the steeper power law, p, 
rather than a decrease in the central flat radius, Rfi. As discussed in Section 4.2, there is relatively little difference in 
the quality of fit (comparing typical values) for the modified isothermal model when a background column density 
is or is not included, however, the inclusion of a background column density terms makes a substantial difference in 
the quality of fits for the purely isothermal model. 

ACKNOWLEDGEMENTS 

We thank the referee for a thoughtful and thorough review which helped to strengthen our paper. MK thanks Philip 
Girichidis for sharing his turbulence generator with us, and Thierry Sousbie for providing us with a pre-release version 
of DisPerSE. HK thanks Rowan Smith for helpful discussions about analyzing filaments in numerical simulations. 

































18 


especially in comparing her recent results with ours. HK thanks Doris Arzoumanian for providing some clarification 
to the radial column density profile fitting method used by the Herschel team. HK also thanks Doug Johnstone and 
Phil Myers for useful discussions about various aspects of this work. HK acknowledges support from the Banting 
Postdoctoral Fellowships program, administered by the Government of Canada. MK acknowledges financial support 
from the National Sciences and Engineering Research Council (NSERC) of Canada through a Postgraduate Scholarship. 
REP is supported by a Discovery Grant from the Natural Sciences and Engineering Research Council (NSERC) of 
Canada. The flash code was in part developed by the DOE NNSA-ASC OASCR Elash Center at the University 
of Chicago. This work was made possible by the facilities of the Shared Hierarchical Academic Research Computing 
Network (SHARCNET: www.sharcnet.ca) and Compute/Calcul Canada. 

REFERENCES 


Alves, J., Lada, C. J., Lada, E. A., Kenyon, S. J., & Phelps, R. 
1998, ApJ, 506, 292 

Andre, P., Men’shchikov, A., Bontemps, S., et al. 2010, A&A, 

518, L102 

Arzoumanian, D., Andre, P., Peretto, N., &: Konyves, V. 2013, 
A&A, 553, A119 

Arzoumanian, D., Andre, P., Didelon, P., et al. 2011, A&A, 529, 
L6 

Bally, J., Lanber, W. D., Stark, A. A., & Wilson, R. W. 1987, 
ApJ, 312, L45 

Banerjee, R., Pudritz, R. E., &: Anderson, D. W. 2006, MNRAS, 
373, 1091 

Banerjee, R., Vazquez-Semadeni, E., Hennebelle, P., & Klessen, 

R. S. 2009, MNRAS, 398, 1082 
Bergin, E. A., & Tafalla, M. 2007, ARA&A, 45, 339 
Beuther, H., Vlemmings, W. H. T., Rao, R., & van der Tak, 

F. F. S. 2010, ApJ, 724, L113 
Boldyrev, S. 2002, ApJ, 569, 841 

Bourke, T. L., Myers, P. C., Caselli, R, et al. 2012, ApJ, 745, 117 
Contreras, Y., Rathborne, J., & Garay, G. 2013, MNRAS, 433, 
251 

Crutcher, R. M., Wandelt, B., Heiles, C., Falgarone, E., Sz 
Troland, T. H. 2010, ApJ, 725, 466 
Falgarone, E., Troland, T. H., Crutcher, R. M., & Paubert, G. 

2008, A&A, 487, 247 

Federrath, C., Banerjee, R., Clark, P. C., Sz Klessen, R. S. 2010, 
ApJ, 713, 269 

Federrath, C., Klessen, R. S., & Schmidt, W. 2008, ApJ, 688, L79 
Fernandez-Lopez, M., Arce, H. G., Looney, L., et al. 2014, ApJ, 
790, L19 

Fiege, J. D., & Pudritz, R. E. 2000, MNRAS, 311, 85 
Fischera, J., & Martin, P. G. 2012, A&A, 542, A77 
Foster, J. B., Rosolowsky, E. W., Kauffmann, J., et al. 2009, ApJ, 
696, 298 

Fryxell, B., Olson, K., Ricker, P., et al. 2000, ApJS, 131, 273 
Girart, J. M., Beltran, M. T., Zhang, Q., Rao, R., &: Estalella, R. 

2009, Science, 324, 1408 

Girichidis, P., Federrath, C., Banerjee, R., Sz Klessen, R. S. 2011, 
MNRAS, 413, 2741 

Gomez, G. C., &: Vazquez-Semadeni, E. 2014, ApJ, 791, 124 
Goodman, A. A., Pineda, J. E., Sz Schnee, S. L. 2009, ApJ, 692, 
91 

Hacar, A., & Tafalla, M. 2011, A&A, 533, A34 
Hacar, A., Tafalla, M., Kauffmann, J., Sz Kovacs, A. 2013, A&A, 
554, A55 

Hatchell, J., Wilson, T., Drabek, E., et al. 2013, MNRAS, 429, 
LIO 

Heitsch, F. 2013a, ApJ, 769, 115 
—. 2013b, ApJ, 776, 62 

Hennebelle, R, & Andre, P. 2013, A&A, 560, A68 
Hennemann, M., Motte, F., Schneider, N., et al. 2012, A&A, 543, 
L3 

Henshaw, J. D., Caselli, P., Fontani, F., Jimenez-Serra, L, Sz Tan, 
J. C. 2014, MNRAS, 440, 2860 
Heyer, M. H., & Brunt, C. M. 2004, ApJ, 615, L45 
Hill, T., Andre, P., Arzoumanian, D., et al. 2012, A&A, 548, L6 
Holland, W. S., Bintley, D., Chapin, E. L., et al. 2013, MNRAS, 
430, 2513 

Inutsuka, S.-L, Sz Miyama, S. M. 1997, ApJ, 480, 681 
Johnstone, D., Di Francesco, J., Sz Kirk, H. 2004, ApJ, 611, L45 
Juvela, M., Malinen, J., Sz Lunttila, T. 2012a, A&A, 544, A141 
Juvela, M., Ristorcelli, L, Pagani, L., et al. 2012b, A&A, 541, A12 
Kauffmann, J., Pillai, T., Shetty, R., Myers, P. C., Sz Goodman, 


Kirk, H., Johnstone, D., Sz Di Francesco, J. 2006, ApJ, 646, 1009 
Kirk, H., Myers, P. C., Bourke, T. L., et al. 2013, ApJ, 766, 115 
Kirk, H., Pineda, J. E., Johnstone, D., Sz, Goodman, A. 2010, 

ApJ, 723, 457 

Klassen, M., Pudritz, R. E., Sz Peters, T. 2012, MNRAS, 421, 

2861 

Klessen, R. S., Heitsch, F., Sz Mac Low, M.-M. 2000, ApJ, 535, 
887 

Konyves , V., Andre, P., Schneider, N., et al. 2013, 

Astronomische Nachrichten, 334, 908 
Konyves, V., Andre, P., Men’shchikov, A., et al. 2010, A&A, 518, 
L106 

Lada, C. J., Alves, J., Sz Lada, E. A. 1999, ApJ, 512, 250 
Larson, R. B. 1981, MNRAS, 194, 809 

MacNeice, P., Olson, K. M., Mobarry, C., de Fainchtein, R., Sz 
Packer, C. 2000, Computer Physics Communications, 126, 330 
Malinen, J., Juvela, M., Rawlings, M. G., et al. 2012, A&A, 544, 
A50 

Markwardt, C. B. 2009, in Astronomical Society of the Pacific 
Conference Series, Vol. 411, Astronomical Data Analysis 
Software and Systems XVHI, ed. D. A. Bohlender, D. Durand, 
Sz P. Dowler, 251 

Moeckel, N., Sz Burkert, A. 2014, ArXiv e-prints 

Motte, F., Zavagno, A., Bontemps, S., et al. 2010, A&A, 518, L77 

Mouschovias, T. C., Sz Spitzer, Jr., L. 1976, ApJ, 210, 326 

Myers, P. C. 2009, ApJ, 700, 1609 

—. 2011, ApJ, 735, 82 

Nutter, D., Kirk, J. M., Stamatellos, D., Sz Ward-Thompson, D. 

2008, MNRAS, 384, 755 

Offner, S. S. R., Klein, R. L, McKee, C. F., Sz Krumholz, M. R. 

2009, ApJ, 703, 131 

Olson, K. M., MacNeice, P., Fryxell, B., et al. 1999, in Bulletin of 
the American Astronomical Society, Vol. 31, American 
Astronomical Society Meeting Abstracts, 1430—h 
Ostriker, J. 1964, ApJ, 140, 1056 

Palmeirim, P., Andre, P., Kirk, J., et al. 2013, A&A, 550, A38 
Peretto, N., Andre, P., Konyves, V., et al. 2012, A&A, 541, A63 
Peretto, N., Fuller, G. A., Andre, R, et al. 2014, A&A, 561, A83 
Pineda, J. E., Goodman, A. A., Arce, H. G., et al. 2011, ApJ, 

739, L2 

Pirogov, L. E. 2009, Astronomy Reports, 53, 1127 
Plummer, H. C. 1911, MNRAS, 71, 460 
Recchi, S., Hacar, A., Sz Palestini, A. 2013, A&A, 558, A27 
—. 2014, MNRAS, 444, 1775 

Rosolowsky, E. W., Pineda, J. E., Foster, J. B., et al. 2008, ApJS, 
175, 509 

Schnee, S., Rosolowsky, E., Foster, J., Enoch, M., Sz Sargent, A. 
2009, ApJ, 691, 1754 

Schneider, N., Csengeri, T., Hennemann, M., et al. 2012, A&A, 
540, Lll 

Schneider, S., Sz Elmegreen, B. G. 1979, ApJS, 41, 87 
Seifried, D., Banerjee, R., Klessen, R. S., Duffin, D., Sz Pudritz, 

R. E. 2011, MNRAS, 417, 1054 

Smith, R. J., Glover, S. C. O., Sz Klessen, R. S. 2014, MNRAS, 
445, 2900 

Sousbie, T. 2011, MNRAS, 414, 350 

Sousbie, T., Pichon, C., Sz Kawahara, H. 2011, MNRAS, 414, 384 
Walawender, J., Bally, J., Francesco, J. D., Jprgensen, J., Sz 
Getman, K. . 2008, NGC 1333: A Nearby Burst of Star 
Formation, ed. B. Reipurth, 346 


